Impact of Cholesterol on Voids in Phospholipid Membranes 
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Free volume pockets or voids are important to many biological processes in cell membranes. Free volume 
fluctuations are a prerequisite for diffusion of lipids and other macromolecules in lipid bilayers. Permeation of 
small solutes across a membrane, as well as diffusion of solutes in the membrane interior are further examples 
of phenomena where voids and their properties play a central role. Cholesterol has been suggested to change the 
structure and function of membranes by altering their free volume properties. We study the effect of cholesterol 
on the properties of voids in dipalmitoylphosphatidylcholine (DPPC) bilayers by means of atomistic molecular 
dynamics simulations. We find that an increasing cholesterol concentration reduces the total amount of free 
volume in a bilayer. The effect of cholesterol on individual voids is most prominent in the region where the 
steroid ring structures of cholesterol molecules are located. Here a growing cholesterol content reduces the 
number of voids, completely removing voids of the size of a cholesterol molecule. The voids also become 
more elongated. The broad orientational distribution of voids observed in pure DPPC is, with a 30 % molar 
concentration of cholesterol, replaced by a distribution where orientation along the bilayer normal is favored. 
Our results suggest that instead of being uniformly distributed to the whole bilayer, these effects are localized 
to the close vicinity of cholesterol molecules. 



I. INTRODUCTION 

In addition to space occupied by atoms of, e. g., lipid, sterol, 
and protein molecules, cell membranes'-^'^ incorporate free 
space ox free volume ^/^/^^^^^^ Membranes may be thought 
of as complex porous structures with a varying number of free 
volume pockets or voids. The voids might assume different 
sizes, shapes, and orientations. The nature of the voids is dy- 
namic: they can be generated or annihilated by trans-gauche 
isomerizations in the hydrocarbon tails of lipids molecules, or 
less frequently, by the movement of whole lipidsi^ 

Voids are crucial for dynamic processes in membranespiS 
Local free volume fluctuations in bilayers are supposed to 
give rise to jump diffusion of the lipids that constitute the 
membrane bilayer "^'^ Another important process where voids 
play a central role is the diffusion of small solutes within 
or across membranes. For instance the electron transport in 
mitochondria and chloroplasts is primarily influenced by the 
availability of voids for the diffusion of the electron carry- 
ing quinone molecule pii Passive transport of water, oxygen, 
small organic molecules, and small ions to and from the cell 
across the plasma membrane is an important process involv- 
ing voids. ^iiSi^ii^iiiii^ii^ Further, general anesthesia might be 
partially explained by changes in the packing and void dis- 



tribution of membranes caused by the partitioning and diffu- 
sion of anesthetics in membranes. These changes have been 
suggested to lead to modifications in the lateral pressure pro- 
file of the membrane, and will consequently alter the structure 
and function of the membrane proteins. '^ '^ Voids are also im- 
portant for the interpretation of fluorescence anisotropy mea- 
surements. Decreased fluorescence anisotropy caused by the 
increased mobility of fluorescent probes such as diphenyl- 
hexatriene (DPH) or pyrene could be due to generation of 
voids. '" '^ Finally, voids are believed to be involved in thin- 
ning transitions in smectic membranesi^ 

There have been several attempts to formally describe the 
effect of voids on the dynamic processes in cell membranes or 
the simpler one- and multi-component lipid bilayers. The first 
free volume and free area theories appeared decades ago4iiSl 
Free area theories in particular have been used for interpreting 
the mobilities of lipid or sterol molecules in model bilayers.iii 
The basis of these theories is that a particle attempting a jump 
needs, in addition to a sufficient activation energy, a large 
enough empty site to jump to. The diffusion coefficient is 
therefore thought to depend on the ratio of the close-packed 
area of the diffusing particle and the average free area per 
particle.'* 

Free area theories have certain shortcomings. They are two- 
dimensional, i. e., it is assumed that the free volume properties 
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do not vary significantly in the direction of the bilayer nor- 
mal. Unfortunately, bilayers are remarkably heterogeneous in 
the normal direction, and the amount of free volume varies 
strongly across a bilayer. ^-^^ More important, free area theo- 
ries are mean-field, that is, all predictions are based on the 
average free area available per lipid. It is impossible to dis- 
tinguish between very few large voids and numerous tiny 
ones. Such a distinction should be highly relevant from the 
point of view of, e.g., solute diffusion; unhindered motion 
within a substantial void is expected to differ from jumps be- 
tween isolated voids. Hence, free area theories should not 
be used for quantitative predictions. Of late, there have been 
a few attempts to develop more sophisticated free volume 
theoriesi22i2i These analytical approaches have enhanced our 
understanding of the impact of voids on dynamic processes. 

A detailed understanding of the microscopic properties of 
voids and their role in dynamical processes will, however, not 
emerge from analytical considerations, but can be acquired 
using computer simulations. In two early molecular dynamics 
(MD) studies, Bassolino-Klimas et al^/^ focused on the diffu- 
sion mechanisms of benzene molecules in phospholipid bilay- 
ers at different temperatures in the fluid phase. The size distri- 
bution of voids was used to explain the motion of the benzene 
molecules. Further, the size distribution was found to depend 
on temperature. Marrink et al."^ performed a detailed anal- 
ysis of the properties of voids in a dipalmitoylphosphatidyl- 
choline (DPPC) bilayer in the fluid phase. The size, shape, 
and orientation of the voids were found to vary significantly 
in the direction of the bilayer normal. More recently, Shinoda 
et al.^* considered the effect of chain branching on the per- 
meability of lipid bilayers, and briefly addressed the issue of 
voids. They concluded that chain branching would seem to 
reduce the probability of observing large voids. Jedlovszky et 
alpi^ have briefly considered the impact of cholesterol on the 
number of voids in dimyristoylphosphatidylcholine (DMPC) 
bilayers. In all, computer simulations can be considered a use- 
ful tool for studying voids. 

In this computational study, we focus on the impact of 
cholesterol on the properties of voids in lipid bilayers. Choles- 
terols are rigid molecules that seem to enhance the proper- 
ties of cell membranes and allow for wider variations in the 
lipid composition of the membrane. Important effects of a 
finite cholesterol concentration on bilayers in the fluid state 
include changes in passive permeability of small solutes 
and suppressed lateral diffusion of phospholipids in bilayers 
containing cholesterol. "^'^ '^^i^^i^^i-^" These effects suggest that 
cholesterol occupies free volume in lipid bilayers. "^'^-^^ The 
void-filling capacity of cholesterol is also thought to be rele- 
vant for the formation of membrane domains enriched in gly- 
cosphingolipids and cholesterol.^' Our intention is to study in 
microscopic detail how the presence of cholesterol molecules 
influences the properties of voids in phospholipid bilayers. 

We have performed 100 ns MD simulations of 
DPPC /cholesterol bilayers in the fluid state, with cholesterol 
concentrations ranging from to 29.7 mol %.^^ The first step 
in characterizing the voids is to single out the free volume by 
mapping the bilayer on a rectangular three-dimensional grid 
and checking which grid points lie outside the van der Waals 



radii of any atoms in the system. The individual free volume 
pockets or voids are identified by using a weighted tree-based 
union /find algorithm with path compression adapted from 
Ref. |32. Finally, the properties of the voids are analyzed 
with the aid of principal component analysis, a technique 
that yields information on the shape and orientation of voids. 
Combining these methods, we investigate the effects of 
cholesterol, permeant size, and location in the bilayer on the 
properties of voids. We first focus on the size distribution of 
voids, considering also percolation of free volume. The shape 
and orientation of voids, and the possible size dependence 
of these properties, are also examined. Finally, we look into 
local effects on voids in the close vicinity of cholesterol 
molecules. 

These studies are qualitative rather than quantitative in na- 
ture. Our intention is to identify trends and describe new phe- 
nomena, simultaneously attempting to shed light on the mi- 
croscopic reasons behind them. 

Our results suggest that replacing phospholipid molecules 
by cholesterols reduces the total volume of the bilayer In 
addition, both the total free volume and the free volume frac- 
tion decrease with an increased cholesterol content. As for the 
void distribution, the effects of cholesterol on the properties of 
voids are most prominent in the region where the rigid steroid 
ring parts of cholesterol molecules reside. Other parts of the 
bilayer are affected to a lesser degree. In the region with the 
ring structures, the presence of cholesterol reduces the num- 
ber of voids, completely removing voids of its own size. The 
voids become more elongated and instead of a broad orien- 
tational distribution, which is observed in pure phospholipid 
bilayers, with 29.7 % cholesterol voids are preferentially ori- 
ented along the bilayer normal. We also find that the effect of 
cholesterol is local, i. e., the effects are most apparent in the 
close vicinity of cholesterol molecules. 



II. MODEL AND SIMULATION DETAILS 

The model systems we studied are fully hydrated lipid bi- 
layers consisting of 128 macromolecules, i. e., DPPCs and 
cholesterols, and 3655 water molecules. In this contribution 
we focus on a pure DPPC bilayer and DPPC /cholesterol bi- 
layers with five different cholesterol molar fractions: x — 
%, 4.7 %, 12.5 %, 20.3 %, and 29.7 %. The underlying MD 
simulations have been described and validated elsewhere,— 
and in the following we shall merely summarize the most im- 
portant details. The previously published results have been 
shown to be in good agreement with experimental studiesi2^ 

The united atom force fields for the DPPC and choles- 
terol molecules were adopted from earlier studiesi^i^^f^ As 
an initial configuration for the pure DPPC bilayer we used 
the final structure of run E in Ref. '3^. For systems with fi- 
nite cholesterol concentrations, the initial configurations were 
constructed by replacing randomly selected DPPC molecules 
with cholesterols such that same number of DPPC molecules 
was replaced in each monolayer To fill the free volume left by 
replacing DPPC molecules by the somewhat smaller choles- 
terols, the system was thoroughly equilibrated, see Ref. |23 
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for details. The bilayer was aligned such that it lies in the xy 
plane. The bilayer normal is thus parallel to the z-axis. 

The MD simulations were run at a temperature T ~ 323 K 
using the GROMACS^^ molecular simulation software. Af- 
ter the initial equilibration, the system was kept in the NpT 
ensemble at a pressure of 1 bar using the Berendsen barostat 
and thermostat.^^ The duration of each simulation was 100 ns 
and the time step was chosen to be 2.0 fs. For long-range elec- 
trostatic interactions we used the Particle-Mesh Ewald (PME) 
method)^ which has been shown to yield reliable results for 
bilayers42il 



III. METHODS 

The details of the technique outlined in this section, along 
with a computer code, will be published elsewhere.l^ 



A. Finding Free Volume: Grid Approach 

To locate the free volume, the bilayers were mapped onto 
three-dimensional rectangular grids, an approach outlined in 
Ref. 9. The grid elements outside the van der Waals radii 
of any atomic groups were designated as free volume and 
those inside as occupied volume. The van der Waals radius 
of a given atomic group was taken to be half of the shortest 
distance where the Lennard-Jones interaction of that atomic 
group with an identical atomic group is zero. The Lennard- 
Jones parameters were extracted from the specification of the 
force field. 

The grid thus constructed stretched, in the normal direc- 
tion, between the points where the density of water starts to 
deviate from its bulk value. The number of grid elements in 
each direction was chosen such that the linear size of an el- 
ement was 0.075 nm on the average. As the area of the bi- 
layer decreases and the thickness increases with an increasing 
cholesterol content, the number of elements had to be set 
separately for each system with a different cholesterol con- 
centration. For instance in the case of pure DPPC, the total 
area of the bilayer is 41.9 nm^ and the distance between the 
points where the density of water no longer corresponds to 
the bulk value is 5.8 nm. With 29.7 % cholesterol, the corre- 
sponding values are 26.9 nm^ and 6.4 nm. In the case of pure 
DPPC, the number of grid elements in the x and y directions 
was taken to be 86. The number of elements in the normal or 
z direction was 78. The corresponding values for the bilayer 
with 29.7 % cholesterol were 74 and 82. 

We also used finer grids to make sure the results were not 
influenced by the resolution. In turns out that when the linear 
size of a grid element is reduced to 0.05 nm, all quantities we 
have computed change only slightly. The changes are clearly 
within the error bars originating from a finite simulation time 
and correlated samples. 

In constructing the grids, in the spirit of Ref. 9, we took 
into account the finite size of small solutes diffusing in bilay- 
ers. Due to the finite size, not all grid elements outside the 
van der Waals radii of the atoms can be accessed by the center 



of mass (CM) of the diffusing solute. Therefore, in addition 
to the empty free volume, i. e., all the free volume outside the 
van der Waals radii of the atoms, we also studied the so-called 
accessible free volume. Accessible free volume is calculated 
by adding the van der Waals radius of the diffusing solute to 
the van der Waals radii of the atoms in the bilayer, and corre- 
sponds to the free volume which is accessible to the CM of the 
diffusing solute. Several solute sizes with radii r ranging be- 
tween and 0.2 nm were studied. For comparison, the effec- 
tive van der Waals radii of bare sodium, chloride, and potas- 
sium ions, as well as those of water and oxygen molecules, 
are between 0. 1 and 0.2 nm. The same applies for example to 
the general anesthetic Xenon. 

The resulting grids were used to compute the average free 
volume fraction as a function of the distance from the bilayer 
center. From these it was easy to extract the total free vol- 
ume in each system. The grids were also the starting point for 
studying the properties of voids for each system. 



B. Identifying and Characterizing Voids 

To characterize the distribution of voids, we need to iden- 
tify all voids, that is, all clusters of free volume, in our grid. 
Computing and characterizing the distribution of contiguous 
clusters of occupied or unoccupied sites is a well-known prob- 
lem, usually solved by a union /find approach. It works as 
follows. The sites, in this case the empty grid elements, must 
be traversed. Whenever we arrive at a new site, a find oper- 
ation is performed, i. e., we check which clusters the nearest 
neighbor sites, if unoccupied, belong to. The union operation 
follows: if there are unoccupied nearest neighbor sites that 
belong to different clusters, those clusters are amalgamated. 
Otherwise, no union of two clusters is needed. 

Modern union /find algorithms are weighted, tree-based, 
and use path compression. Each cluster is stored as a tree 
such that the nodes of the tree are the sites of the cluster Ev- 
ery cluster has a root that corresponds to the root node of the 
tree, and the other nodes have pointers to either the root or to 
another node in the same tree. By following a succession of 
such pointers we can locate the root of any tree. The find oper- 
ation therefore consists of traversing trees by following point- 
ers to the root nodes. If the nearest neighbor sites of a new 
site lead to the same root node, the sites belong to the same 
cluster. The find operation can be implemented such that after 
the traversal is completed, the pointers along the path all point 
directly to the root node. This is called path-compression, and 
it speeds up future traversals. The union of two clusters is 
achieved by adding a pointer from the root node of one tree to 
the root of the other Hence, one of the trees becomes a sub- 
tree of the other. This is done in weighted fashion such that 
the smaller tree becomes a subtree of the larger tree, which is 
easily accomplished by storing the size of each cluster in the 
corresponding root node. 

Even though the above algorithm may seem a little com- 
plicated, it can be implemented, using recursion, in a very 
elegant and concise manner Moreover, its performance is in 
most cases far superior-'^ to simple relabeling algorithms such 
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as that introduced in Ref . '44' We chose therefore to adapt the 
weighted tree-based union-find code with path compression 
given in Ref.l32l 

The end resuh of the union / find algorithm for a given con- 
figuration is a tree structure containing the information on how 
the unoccupied sites with known coordinates are connected 
to each other From the tree structure we constructed arrays 
where each array element corresponds to a free volume cluster 
or void and contains the relevant properties of the void. Such 
properties are the volume of the void, the location of its CM, 
and its maximum and minimum coordinates in each direction. 
Further, we stored the covariance matrix extracted from the 
coordinates of the grid elements that belong to the void. The 
covariance matrix can be used to characterize the shape and 
orientation of the void, as described in the next section. 



spherical clusters. For the elongated voids, we monitored the 
angle cf) between the longest axis and the bilayer normal. 

As for the practical details of the PCA calculation, we com- 
puted the eigenvectors and eigenvalues using singular value 
decomposition. Further, we only used PCA for voids that 
were larger than 4 x 10^"^ nm"^ and smaller than 0.13 nm'^. The 
very small voids were rejected, because it is somewhat artifi- 
cial to discuss the shape of an object consisting of one or two 
rectangular grid elements. The large voids were not consid- 
ered, because their shapes were, in most cases, not ellipsoidal, 
but rather more complicated with branches and compartments. 
Applying PCA to these voids would have been meaningless. 

D. Four Region Model 



C. Principal Component Analysis 

The shape and orientation of the voids were studied using 
principal component analysis (PCA).^^ PCA is a versatile sta- 
tistical method, which is widely used in analysis and com- 
pression of multivariate data, signal processing, and neural 
computing. Here we shall confine ourselves to the case of 
three-dimensional data. 

In applying PCA, we view each void as a three-dimensional 
ellipsoidal cloud. The principal component vectors are or- 
thogonal and represent the directions with maximum variabil- 
ity. In other words, the principal components define the axes 
of the ellipsoidal void such that the first principal component 
is the longest axis, the second is the next longest axis, and the 
third is the shortest axis. 

The principal components may be computed as follows. 
For a void consisting of N grid elements with coordinates 
{rijfLi, fii's'^ compute the mean (or CM) R and the co- 
variance matrix C. The diagonal of the covariance matrix 
contains the variances of the Cartesian components of r^. The 
off-diagonal matrix elements represent correlations between 
the different components of r^. The covariance matrix is sym- 
metric, and by finding its eigenvalues and eigenvectors we find 
an orthogonal set of basis vectors. By ordering the eigenval- 
ues, we obtain an orthogonal basis with the first eigenvector 
having the direction of the largest variance of the data and so 
on. This is the set of uncorrected principal components, and 
simultaneously, the axes of our ellipsoidal void. The square 
roots of the eigenvalues, in turn, are proportional to the lengths 
of the axes. 

The shape of a void can be deduced from the eigenvalues 
af of the covariance matrix. We chose to focus on the ratios 
of the square roots of the eigenvalues, denoted by ai / <Tj . In 
case iTi/ct2 ~ 0-2 /c3 ~ 1, the void is spherical, otherwise it is 
elongated. The principal components were used to calculate 
the orientation of each elongated or non-spherical void. We 
chose to regard all voids with ai/a^ < 1.5 as spherical. Even 
though this is not a very strict definition of sphericity, most 
voids appeared to be non-spherical. Further, this definition 
was found to be robust, since changing the somewhat arbitrary 
1.5 to lower values had virtually no effect on the number of 



The bilayer is highly heterogeneous in the normal direc- 
tion, and both the free volume fraction and the properties of 
the voids are anticipated to vary with the position along the bi- 
layer normal.«2.It is therefore not sensible to study the void dis- 
tribution averaged over the whole bilayer, but to discretize the 
bilayer into regions with more homogeneous compositions. 
The properties of voids can then be studied in these separate 
regions. We chose to analyze the data in terms of the four 
region model introduced by Marrink et al.^ '^ 

The original partition into four regions was slightly mod- 
ified to better probe the effects of cholesterol on the proper- 
ties of voids, the main modification being that Region 3 was 
taken to be the part of the bilayer where the density of choles- 
terol steroid ring structures is high. We expect the effects of 
cholesterol on the voids to be most prominent in this region, 
as the steroid rings rigidify the disordered tails of the DPPC 
molecules, and enhance the molecular packingi2Si2£ 
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FIG. 1: Mass density profiles of molecules and atomic groups at 
X ~ 20.3 %. The mass density of the cholesterol hydroxyl groups 
has been scaled by a factor of ten. The division into four regions, 
indicated here by numbers 1—4, is based on the mass density profiles, 
see Sect. lIIIDI 

Figure \l\ illustrates how we have divided the bilayer into 
four distinct regions. Region 1 ranges from the point where 
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the density of water starts to deviate from the bulk value to the 
point where the densities of DPPC and water are equal. This 
region therefore contains mostly water molecules and a low 
density of DPPC headgroups. Region 2 extends from the point 
where the densities of water and DPPC coincide to the point 
where the density of cholesterol ring structures is half of its 
maximum plateau value. This region contains mainly parts 
of DPPC molecules: the choline, phosphate, and glycerol den- 
sities peak in this region, and there is also a finite density of 
DPPC tails. Further, the hydroxyl groups of cholesterol reside 
mainly in this region. Region 3 is defined between the points 
where the cholesterol steroid ring density assumes half of its 
maximum value, and contains in addition the middle parts of 
the DPPC tails. The remaining part of the bilayer is Region 4, 
i. e., this region is in the center of the bilayer. Here we expect 
to find a low density of DPPC and cholesterol tails, with the 
end methyl groups being mainly located in this region. 

The boundaries between the regions vary with cholesterol 
concentration. This is due to the fact that an increasing choles- 
terol content changes the densities of the various components 
such that they are shifted away from the bilayer center.^^ 
Hence, in the case of pure DPPC, the boundaries from the 
outer boundary of Region 1 to the boundary between Regions 
3 and 4 were set to 2.9 nm, 2.1 nm, 1.4nm, and 0.6 nm from 
the bilayer center With 29.7 % cholesterol, the corresponding 
values were 3.2 nm, 2.5 nm, 1.7 nm, and 0.8 nm. 



E. On Percolation Theory 

Percolation theory can be used to predict certain proper- 
ties of clusters of occupied or unoccupied sites or bonds in 
a lattice.'*^ The applications of percolation theory range from 
modeling forest fires or the distribution of oil and natural gas 
in the porous rock in oil reservoirs to studying, say, voids in 
lipid bilayers. The properties of the clusters are studied as 
a function of the fraction of occupied or unoccupied sites or 
bonds, depending on the application. In our case, we are in- 
terested in the free volume fraction 
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{Vtc 



(1) 



where (Vfioc) is the average free volume in a bilayer and (Vtot) 
is the total volume of the bilayer An important free volume 
fraction is the percolation threshold pc, at which there is a free 
volume cluster stretching in the x, y, or z direction from one 
side of the system to the opposite side. The size distribution 
of clusters depends heavily on p and the vicinity of the per- 
colation threshold. At p > pc, there is a free volume cluster 
that spans the whole bilayer, and a large fraction of all empty 
grid elements belong to this percolating cluster. When p < pc, 
there is no such percolating cluster. 

The detailed form of the size distribution of clusters can- 
not, in most cases, be predicted analytically.'*^ A general scal- 
ing form for the size distribution has been postulated, and the 
non-universal details such as the form of the scaling function 
can be extracted from computer simulations. In certain lim- 



iting cases, however, the exact form of the scaling function 
becomes irrelevantii^ 

Very close to the percolation threshold pc or for very small 
clusters, the size distribution N{V) scales with the cluster vol- 



ume as 



N{v) - y-^. 



(2) 



Here r is a universal scaling exponent, which in three dimen- 
sions assumes the value 2.18."*^ This scaling law is valid for 
clusters whose linear size is smaller than the correlation length 
^, which can be thought of as proportional to a typical cluster 
diameter. At the percolation threshold ^ approaches infinity, 
and Eq. (|2j is valid for all cluster sizes. 

Away from pc and for large clusters the size distribution 
scales as follows:'*'' 



N{V) 



V'^exp{-c"V) P<Pc, 
V~'^' exp(-c"'V^i-i/'') p > Pc. 



(3) 



Here d is the dimensionaUty of the system, c" and cJ" are 
proportional to p — pc, and 9 and 9' are universal exponents. 
In three dimensions 9 — 3/2 and 9' = —1/9.'*^ This scaling 
form applies for clusters whose linear size is larger than ^, in 
practice for almost all clusters in systems that are not in the 
near vicinity of pc- Because of the p — Pc dependence of c" 
and c'", further away from pc the crossover from power-law 
to exponential behavior is observed at smaller values of V. 

We have computed size distributions for systems with dif- 
ferent cholesterol concentrations and penetrant sizes. In in- 
terpreting the size distribution data, it is useful to know how 
close to the percolation threshold a given system is. As op- 
posed to infinite systems, in a finite system where the linear 
dimension in the z direction cannot be increased, only appar- 
ent percolation thresholds can be calculated.^ Computing ap- 
parent percolation thresholds is not a problem, since this is 
the threshold a permeant solute will experience.^ The apparent 
percolation threshold can be defined as the free volume frac- 
tion at which at least 50 % of all systems percolate.^ Alterna- 
tively, Pc can be extracted from maxima in the second moment 
of the cluster size distribution, computed over all clusters but 
the percolating one."*^ 

Unfortunately, in bilayer systems we cannot automatically 
expect Eqs. (|2j and (|3j with three-dimensional critical expo- 
nents to describe the behavior of the size distribution at and 
off the apparent percolation threshold. Predicting the form 
of the size distribution is in fact much more difficult for bi- 
layer systems than for isotropic bulk systems in two or three 
dimensions. First, the bilayer is heterogeneous in the normal 
direction, and consequently, the free volume fraction varies 
with the position along the bilayer normal, i. e., p — p{z), see 
Fig. 13 The bilayer center or Region 4, with the highest free 
volume fraction in the whole bilayer, ^-^^ could be well above 
the percolation threshold, while the other regions are well be- 
low Pc- This problem can be in part remedied by consider- 
ing the four regions with approximately constant free volume 
fractions. Secondly, the z direction is very different from the x 
and y directions, and there might well be different thresholds 
for percolation in different directions. Finally, there are finite 
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FIG. 2: Bilayer volume as function of cholesterol concentration. As 
there is no unique definition for the thickness of a bilayer, the thick- 
nesses have been measured in three ways: between the points where 
the mass density of water starts to deviate from the bulk value (□), 
between the points where the mass density of phospholipids and wa- 
ter are equal (o), and between the maxima in the total electron densi- 
ties (•). 



size effects. The dimension in the normal direction will sel- 
dom exceed 10 nm in bilayer systems. This makes the system 
quasi-two-dimensional rather than three-dimensional. 



IV. RESULTS 
A. Total Volume and Free Volume 

It is believed that the total free volume in a lipid bilayer 
decreases when the cholesterol concentration increases. "^'^ On 
the other hand, it is often assumed that the total volume of the 
bilayer is kept approximately constant with a growing choles- 
terol concentration. This is attributed to the decreasing area 
and increasing thickness of the bilayer. As a cholesterol occu- 
pies less volume than a DPPC, the volume occupied by phos- 
pholipids and cholesterols must diminish with an increasing 
cholesterol concentration. Water is pushed away from the bi- 
layer center,^^ and hence we expect that the volume occupied 
by water molecules should at least not grow. The constant 
total volume and a reduced volume occupied by cholesterol, 
DPPC molecules, and water would together suggest that the 
free volume should, contrary to the general belief, increase. 
To determine whether the total free volume decreases or in- 
creases, we have computed the total volume and total free 
volume as functions of cholesterol concentration. 

The total volume is computed by multiplying the average 
area of the bilayer by the thickness of the bilayer. Because of 
the rough interface between the bilayer and the water phase, 
there is no unique definition for the thickness of the bilayer. 
As measures of thickness we have used the distance between 
the points where the mass density of water starts to deviate 
from the bulk density, between the points where the mass den- 



sities of phospholipids and water coincide, and between the 
maxima in the total electron densities. The total volume as a 
function of the cholesterol concentration is shown in Fig.|2] It 
is clear that regardless of how we define the thickness of the 
bilayer, the total volume decreases with a growing cholesterol 
concentration. 

Unfortunately, there is little experimental evidence to ei- 
ther confirm or contradict this result. To our knowledge, there 
are no systematic studies on the effect of cholesterol on the 
phospholipid bilayer area or volume. Monolayer studies do 
exist/^ but the applicability of these results to bilayers is 
questionable. 

As the total volume and the occupied volume both decrease 
with an increasing cholesterol concentration, it is not immedi- 
ately clear how the free volume should be expected to behave. 
This issue can be resolved by extracting the total free volumes 
from the three-dimensional grids introduced in Sect. IIII Al 
The procedure can be thought of as integrating over the to- 
tal free volumes in thin slices in the xy plane. The range of 
integration must be the thickness of the bilayer, defined as in 
the case of total volume. The results are shown in Fig. 13 It 
is clear that the total free volume decreases with an increasing 
cholesterol concentration, notwithstanding how the thickness 
is measured. We can therefore conclude that both the occupied 
and free volume decrease, together resulting in a reduction of 
the total volume of the bilayer. 
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FIG. 3: Empty free volume as function of cholesterol concentration. 
The bilayer thicknesses have been defined as the distance between 
the points where the density of water starts to deviate from the bulk 
value (□), the distance between the points where the density of phos- 
pholipids and water are equal (o), and the distance between the max- 
ima in the total electron densities (•). The errors are of the order of 
a few per cent. 

Experimental findings are in favor of the reduction of free 
volume. Galla et al.^ have reached this conclusion by com- 
paring their phospholipid diffusion coefficients to the predic- 
tions of free area theory. Almeida et al.'* present empirical 
laws for the behavior of specific volumes and thicknesses 
of DMPC / cholesterol bilayers as functions of temperature. 
Based on these laws they have computed the average areas 
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per molecule at different cholesterol concentrations. Com- 
bining these with close-packed areas for DMPC and choles- 
terol molecules, they have estimated the average free areas per 
phospholipid. The free areas thus obtained explain nicely the 
diffusion coefficients from fluorescence recovery after photo- 
bleaching (FRAP) experimentsii 



B. Free Volume Fraction Profile 
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FIG. 4: Free volume fractions as functions of position along bi- 
layer normal, (a) Empty free volume fraction and accessible free 
volume fractions for bilayer with 20.3 % cholesterol. The penetrant 
radii from top to bottom are 0.00 nm, 0.05 nm, 0.07 nm, 0. 10 nm, and 

0. 15 nm. (b) Accessible free volume fraction for different cholesterol 
concentrations for a penetrant with a radius r = 0.05 nm. In both 
panels the value z — corresponds to the bilayer center, and the 
errors are of the order of a percent. 

Having concluded that cholesterol reduces the total free 
volume, let us now look at the behavior of free volume in 
different parts of the bilayer Free volume fraction profiles, 

1. e., free volume fractions as functions of the position z along 
the bilayer normal, are displayed in Fig. 0] Part (a) of the 
figure illustrates the effect of penetrant size on the free vol- 
ume fraction with 20.3 % cholesterol. The uppermost curve 
is the empty free volume fraction, while the bottom one cor- 
responds to the accessible free volume fraction for penetrants 
with a radius r — 0.15 nm. Part (b) shows the effect of choles- 



terol concentration on the accessible free volume fraction with 
r — 0.05 nm. 

All free volume fraction profiles in Fig.|4] regardless of the 
penetrant size or cholesterol concentration, have maxima in 
the bilayer center and minima in Region 2. The penetrant size 
has a strong effect on free volume fraction: at least a third 
of the total volume at each z seems to be empty free volume, 
while for penetrants with r ~ 0.15 nm, apart from the bi- 
layer center, there is virtually no accessible free volume at 
all. These findings are in good qualitative agreement with the 
results of Marrink et al.,'' considering that those simulations 
were performed at a temperature T = 350 K and that tem- 
perature has a significant effect on free volume.^ Cholesterol 
influences the shape of the profile, reducing the free volume 
fractions in Regions 1 and 3. The decrease in Region 3 must 
be due to the presence of the rigid steroid ring structure. In 
Region 1 , the origin of the effect may be the reorientation of 
the DPPC headgroupsi^ 



C. Percolation 
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FIG. 5: Second moments of void size distributions computed over 
all but percolating cluster as functions of penetrant radii. The open 
circles are for x = % and the filled ones for x ~ 29.7%. The 
errors are of the order of a few percent. 

Examining the nature and location of the percolation tran- 
sition is essential for a proper interpretation of the size distri- 
butions. We have located the percolation threshold by com- 
puting the second moment of the size distribution of voids, 
see Sect. llIIEl The second moments for x — 0% and 29.7 % 
as functions of penetrant radii are shown in Fig. [S] Since a 
maximum indicates percolation, we can identify two percola- 
tion transitions for each cholesterol concentration. For small r 
percolation takes place in x, y, and z directions. In between the 
two maxima the system percolates in the plane of the bilayer, 
but no longer in the normal direction. This xy percolation al- 
ways takes place in Region 4. Both of these transitions are 
important from a biological point of view: the xyz percola- 
tion is relevant for permeation across the bilayer, whereas xy 
percolation in the bilayer center could be significant, e. g., for 
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diffusion of quinone within the bilayer. Finally, at large values 
of r there is no percolation. 

The location of the transitions depends weakly on choles- 
terol concentration. The xyz percolation takes place at r = 
0.048(2) nm when x = 0% and at r = 0.044(2) nm when 
X — 29.7%. The corresponding values for the xy percolation 
are r = 0.080(2) nm and r — 0.073(2) nm. A larger fraction 
of cholesterol therefore seems to imply that slightly smaller 
penetrants perceive the system as percolating. 

As for percolation thresholds, it is meaningless to define 
anything such for the xyz percolation, since the free volume 
fraction varies strongly along the bilayer normal. The xy per- 
colation in the more homogeneous Region 4 is quite another 
matter. In this case we can, using plots similar to Fig.|3 map 
the penetrant sizes for which percolation takes place to free 
volume fractions. The threshold for xy percolation is p = 0.06 
at all cholesterol concentrations. Cholesterol does not alter 
the percolation threshold itself, but it changes slightly the free 
volume fraction for a given penetrant size. This causes the 
weak X dependence of the location of the xy percolation seen 
in Fig. |5l The value of 0.06 is in good agreement with the 
corresponding threshold calculated by Marrink et al. for pure 
DPPC and thus equals their value for soft polymer as welli^ 



D. Size Distributions 

The size distributions are unnormalized to better illustrate 
the behavior of the absolute number of voids in a bilayer of 
a given composition. Therefore N{V) is the average number 
of voids of volume V. As the four regions have very distinct 
material properties,^' '"^ '^ we have considered the size distri- 
butions separately in each region. Further, the vicinity of the 
percolation transitions is significant for the form of the distri- 
bution. We have identified two cases that represent the most 
biologically relevant regimes: r — 0.05 nm at which xy per- 
colation takes place in Region 4 and r = 0.09 nm with no 
percolation. The regime where xyz percolation occurs is quite 
similar to the case of r = 0.05 nm. Besides, even though the 
results presented here are qualitative rather than quantitative 
in nature, describing the qualitative effect of cholesterol in- 
stead of claiming to establish accurate numerical values, pen- 
etrant radii smaller than 0.05 nm are unphysical. 

The effect of cholesterol on the void size distributions 
is largely restricted to those parts where the cholesterol 
molecules are located. In Region 1, regardless of r, choles- 
terol has no impact on the size distributions. This is not sur- 
prising, since no part of cholesterol is present in the water 
phase. The effects are still very minor in Region 2, although 
the hydroxyl groups of cholesterol are mainly located in this 
region. The situation in Regions 3 and 4, where the choles- 
terol steroid ring structures and tails are situated, is illustrated 
in Fig. |6] These are somewhat coupled, since a particularly 
large cluster could extend to both regions, but is assigned to 
Region 4, as its CM is located there. 

The behavior of these two regions with an increasing 
cholesterol concentration depends on r. The case of r = 
0.05 nm, representing all r below the xyz percolation thresh- 
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FIG. 6: Effect of cholesterol concentration on void size distribution. 

(a) Region 3 and r = 0.05 nm. (b) Region 4 and r — 0.05 nm. 
(c) Region 3 and r — 0.09 nm. (d) Region 4 and r = 0.09 nm. 
Percolation in the plane of the bilayer is seen in panel (b) as a separate 
cluster of points in the right-hand comer. The finite statistics lead to 
relative errors that grow with V. For V < 0.01 nm"^ the errors in 
N{V) are smaller than a percent, and in the range 0.01 imi^ < V < 
0.1 nm'^ smaller than ten percent. If V is of the order of 1 nm'', the 
relative errors may be 100 — 300 %. As the data are shown on a 
loglog scale, this is hardly a problem. 

old and above the xy threshold, is portrayed in panels (a) and 

(b) of Fig.|6l The most notable effect of replacing 29.7 % of 
the DPPC molecules by cholesterols is the decreasing connec- 
tivity of free volume. This is manifested in the size of the per- 
colating cluster, which decreases by a factor of three to V ^ 
5 nm'^, approximately ten times the average close-packed vol- 
ume of a cholesterol molecule, Vchoi ~ 0.459(2) nm^^.^^ 
When connectivity is reduced, the number of non-percolating 
voids of all sizes in Region 4 increases: although the total 
free volume decreases, the number of voids increases approx- 
imately by a factor of three. The effects in Region 3 are more 
minor Here the number of voids from V with radii of 0.1 nm 
to V of the order of a tenth of Vchoi increases very slightly. 
The number of large voids with V close to Khoi decreases 
somewhat. 

The situation is quite different when no percolating cluster 
is present, as illustrated by results computed with r = 0.09 nm 
in panels (c) and (d) of Fig. |6l The effects of cholesterol 
are very pronounced in Region 3. Comparing the case of 
29.7 % cholesterol to pure DPPC we note that there are fewer 
voids of all sizes. The number of voids is reduced at least 
by a factor of three, but especially for larger voids the re- 
duction can be as large as a factor of twenty. The largest 
voids with 0.2 nm"^ ^ V ^ 0.7 nm'^, i.e., voids with sizes 
close to Vchoi, are removed completely. The effects are sim- 
ilar, although more minor, in Region 4. Voids with radii up 
to 0.1 nm are unaffected, while the number of intermediate 
and large voids decreases. Here cholesterol removes the voids 
with 1 nm-^ ^ ^ ^ 3 nm'^, larger than its own close-packed 
size. 

Let us finally analyze Fig. |6] from the point of view of per- 
colation theory. Below the xyz percolation threshold, i. e., for 
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both r = 0.05 nm and r = 0.09 nm, there is no percolating 
cluster in Region 3. If we assume that bulk three-dimensional 
percolation theory holds, in Region 3 we expect to observe 
power-law behavior with an exponent 9 = 1.5 at small and in- 
termediate V and cross-over to exponential behavior at large 
V. This indeed seems to be the case, see panels (a) and (c). 
We also note that the cross-over to exponential behavior oc- 
curs at smaller V for r — 0.09 nm. This is to be anticipated, 
since the larger the penetrant radius, the further away from 
the xyz percolation threshold we are. In the case of Region 4 
with r = 0.09 nm, see panel (d), we are very close to the xy 
percolation threshold. At very small V, power-law behavior 
with an exponent close to 2.1 is observed. The linear sizes of 
these small voids must hence be smaller than the correlation 
length ^. At intermediate and large V, the power-law behav- 
ior can still be observed, with an exponent close to 1.5. Pure 
DPPC, which is closer to the percolation threshold than bilay- 
ers with finite cholesterol concentrations, shows no cross-over 
to exponential behavior, while with finite x we observe expo- 
nential behavior for the largest V. The behavior in Region 4 
above the xy percolation threshold, see panel (b), is difficult 
to interpret, as we are in between two percolation thresholds. 
At small and intermediate V we again observe power-law be- 
havior with an exponent close to 2.2. At larger values of V, 
the situation is less clear. 

Having already mentioned the impact of penetrant size and 
region on the void size distributions, let us discuss them at 
some more length. Void size distributions at x = 20.3 % in 
Region 3 for penetrants with radii ranging from 0.05 nm to 
0.15 nm are displayed in Fig.0^a). We can immediately draw 
the obvious conclusion that large penetrants see fewer and 
smaller voids than small ones. As Region 3 with r > 0.05 nm 
is below any percolation threshold, all cases look very similar: 
power-law behavior with an exponent 1.5 at small and inter- 
mediate V and a cross-over to exponential behavior at larger 
V. The cross-over moves to smaller V with an increasing r, as 
we are getting further away from the xyz percolation thresh- 
old. 

The effect of region on the case of x = 20.3 % and 
r — 0.09 nm is shown in Fig. 0b). The smallest number 
of voids, as well as the smallest voids, can be found in the 
tightly packed headgroup region or Region 2. As expected. 
Region 4 features much larger voids than any other region. 
All curves shown in this figure agree well with the predic- 
tions of the bulk three-dimensional percolation theory below 
the percolation threshold. 



E. Shape 

Using PCA we can extract cti, a2, and 0-3, which are pro- 
portional to the lengths of the principal axes of an ellipsoidal 
void. From these we can compute P(cri/CT2, o'2/(T3), a dis- 
tribution describing the probability that a void has given val- 
ues of (T1/CT2 and (72 /f 3. The distribution has been normal- 
ized such that integration over it gives unity. As visualizing 
and comparing different P{ai/a2, oijcrz) can be difficult, we 
have also computed P{ai/aj) describing the probability of 
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FIG. 7: Effect of penetrant size and region on void size distribution, 
(a) Size distributions for different penetrant sizes for the case x ~ 
20.3 % and Region 3. (b) Size distributions in different regions for 
X ~ 20.3 % and r = 0.09 nm. The relative errors are similar to those 
in Fig.|6| 



finding a void with a given value of (Ti/a-j irrespective of ak, 
k ^ These distributions have been scaled such that inte- 
gration over them gives unity. 

Interpretation of the probability distributions will give us 
information on the shape of the voids for which 4 x 10^"^ nm'^ 
<V < O.lSnm^, see Sect. lmcl If ai/a2 ^ a^joz w 1, the 
void is spherical. In case cti 3> (72 « cts, we are dealing with 
a roughly cigar-shaped void. Finally, if cri « (T2 S> (73, our 
void is essentially disk-like. 

The effect of cholesterol on Piaija^^ 02! '^■i) in Region 3 
with r — 0.05 nm is portrayed in Fig.|8] Cholesterol does not 
influence the shape of the voids in any other region, and other 
penetrant sizes give very similar results (data not shown). Pan- 
els (a) and (b) show the probability distribution for pure DPPC 
and panels (c) and (d) for a bilayer with 29.7 % cholesterol. In 
both cases elongated voids dominate the distributions: spheri- 
cal or nearly spherical voids are rare. The presence of choles- 
terol makes the voids more elongated, i. e., larger values of 
(Ti / (T2 and (T2 jaz occur with a higher probability. 

The elongating effect of cholesterol can be further studied 
by focusing on P((Ti/(T2) and P {(72/(13) shown in panels (a) 
and (b) of Fig. |9] These figures show a systematic average 
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FIG. 8: P(cri/cr2, 0-2/0-3) at X = 0% and X = 29.7% in Region 
3. The penetrant size is 0.05 nm. (a) Surface plot at x = 0.0 %. 
(b) Contour plot at x = 0.0%. (c) Surface plot at x = 29.7%. 
(d) Contour plot at x = 29.7 %. The relative errors are smaller than 
ten percent. 
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FIG. 9: P(o-i/o-j) at different cholesterol concentrations and in 
different regions, (a) P(oi/o-2) in Region 3 for r = 0.05 nm. 
(b) P{(J2/o--j,) in Region 3 for r — 0.05 nm. (c) P{ai/a2) for 
X = 20.3% and r = 0.09 nm. (d) P(o-2/o-3) for x = 20.3% and 
r — 0.09 nm. The relative errors are smaller than two per cent when 
Ci/aj < 3, and between two and five percent for 3 < ai/aj < 6. 



elongation with an increasing cholesterol content: the larger 
the concentration, the more elongated voids we see. 

Figure |9l also depicts the behavior of P((Ti/(T2) and 
P(o'2/ct3) in different regions, see panels (c) and (d). In this 
case the bilayer contains 20.3 % cholesterol and r = 0.09 nm, 
which represents well all penetrant sizes. Regions 1 and 2 are 
very similar to each other Most voids larger than the lower 
cutoff 1^ = 4x10^"^ nm'^ have been taken into account in 
PCA, since the vast majority of voids is smaller than the up- 
per cutoff V — O.lSnm'^, see Fig.0 Region 4 appears to 
differ from Regions 1 and 2 only slightly. This is true for 
small voids with V < O.lSnm^^. We should, however, keep in 



mind that Region 4 also contains larger voids, see Fig.El The 
shapes of these holes are not ellipsoidal but more complex. 
The voids in Region 3 are more elongated than elsewhere. 
When the cholesterol concentration increases, this effect be- 
comes slightly more pronounced. If x = %, all regions look 
approximately similar. 



F. Orientation 

The orientations of the principal axes of the ellipsoidal 
voids can be extracted from PCA. We have measured the co- 
sine of the angle (j) between the longest axis of an elongated 
ellipsoid and the bilayer normal. When the longest axis is ori- 
ented along the bilayer normal (/) = 0. In principle, G [0, tt), 
but due to the symmetry properties of ellipsoids the angles 
(j> and TT — are equivalent. We will therefore confine our- 
selves to G [0,7r/2]. To describe the orientation of the 
voids, we have computed the quantity P(cos0), the proba- 
bility that a given void is oriented such that the cosine of the 
angle between its longest axis and the bilayer normal is cos (j>. 
The probability distribution has been normalized to give unity 
when integrated over, d (cos (/))P (cos (/)) = 1. Should we 
wish to think in terms of cj) rather than cos (j>, averages are com- 
puted by integrating over cj) and the weight of a given angle is 
P(cos (h) sin 6. 
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FIG. 10: Orientational distributions at different cholesterol concen- 
trations and in different regions, (a) P(cos 0) sin (ji in Region 3 
for r — 0.05 nm. (b) P(cos<^) in Region 3 for r = 0.05 nm . 
(c) P(cos (j)) sin </> for X = 20.3 % and r = 0.09 nm. (d) P(cos 0)0 
for X = 20.3 % and r — 0.09 nm. The relative errors are of the order 
of five percent for both P(cos 0) and P(cos 0) sin 0. 

The impact of cholesterol concentration on P(cos0) and 
P(cos 0) sin depends on which region of the bilayer we 
study. There are hardly any changes in either Region 1 or 2. 
The effect of cholesterol is strongest in Region 3, and largely 
independent of r. The case of r = 0.05 nm is portrayed in 
panels (a) and (b) of Fig. ^] In pure DPPC, angles (p > O.Itt 
are favored, while orientation along the bilayer normal hardly 
occurs at all. At x = 29.7 %, the angle with the largest weight 
is = O.OGtt. Orientation along the bilayer normal or close to 
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it is favored, and perpendicular orientation is rare. At choles- 
terol concentrations between and 29.7 %, the orientational 
distributions change systematically from one extreme to the 
other The presence of cholesterol can therefore be said to 
promote orientation along the bilayer normal in Region 3. 

In Region 4 the effects of cholesterol are less pronounced 
and depend on r (data not shown). Between the xy and xyz 
percolation thresholds, e. g., in the case of r = 0.05 nm, ori- 
entation along the bilayer normal is favored. An increasing 
cholesterol concentration seems to accentuate this effect. Be- 
low the xy percolation threshold, for r = 0.09 nm and larger, 
perpendicular orientation is favored and few voids are oriented 
along the bilayer normal, see panel (c) of Fig. [TO] Cholesterol 
has a minor effect here, suppressing the weight of intermedi- 
ate angles cj) « 0.257r and making perpendicular orientational 
slightly more probable. As discussed in the end of Sect. llVFl 
it should be kept in mind that these distributions have been 
computed for voids smaller than V — 0.13 nm"^. For instance 
in the case of r = 0.05 nm. Region 4 contains a large perco- 
lating cluster of a complex shape, oriented more or less in the 
plane of the bilayer 

The situation in different regions at x = 20.3 % and for 
r = 0.09 nm is portrayed in panels (c) and (d) of Fig. 1101 
The effect of changing x and t has been described above. In 
Regions 1 and 2, the orientational distribution of voids is very 
broad with a peak close to </> = 7r/2. Angles near (1)^0 
hardly occur at all. In Region 3, orientation along the bilayer 
normal or close to it is favored. Region 4, as described above, 
contains mostly voids that are oriented perpendicular to the 
bilayer normal. 



G. Size, Shape, and Orientation 

So far we have considered the shape and orientation of all 
voids irrespective of their size. In this section we shall investi- 
gate whether voids of different sizes differ in average shape or 
orientation. For this purpose we have computed P{V, (Ji/(T2), 
the probability of finding a void of volume V with a given 
value of (Ti/fT2, and P(y,cos0), the probability of finding a 
void of volume V with a given cos (j). Both probability distri- 
butions have been normalized such that integration over them 
results in unity. 

Unfortunately, P{V,ai/(T2) and P{V,cos(l)) are quite 
noisy and extremely difficult to visualize. Therefore we have 
calculated the average values of (Ti / (72 and cos </> for each V, 
defined as follows: 
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These might still be somewhat noisy, as the number of voids 
of a given size V may be quite small. 

In analyzing {<7i/a2)v and (cos we shall focus on the 
effects of cholesterol in Region 3 and for r = 0.05 nm. Based 
on our studies of P{<Ji/a2), P(cos cj)), and related quantities. 



we expect the impact of cholesterol to be most significant in 
this region. Further, we do not anticipate our results to change 
significantly with r. 
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FIG. 11: Averages {(Ji/(J2)v as functions of V in Region 3 and for 
r — 0.05 nm. The relative errors grow with V: for V < 0.05 nm'^ 
the errors are smaller than three percent, for 0.05 nm'^ < < 
O.lOnm'^ they grow linearly in the range 4 — 10%, and for larger 
V they can be up to 20 %. 

The behavior of {iJi/a2)v with V is shown in Fig.^] In 
the case of pure DPPC, our data indicate that voids of all 
sizes should have approximately similar shapes. With 29.7% 
cholesterol small voids with V < 14hoi/10 seem to be more 
elongated than larger ones. The shapes of large voids are little 
affected by an increasing cholesterol concentration. Based on 
these data, we can conclude that the presence and amount of 
cholesterol appear to mainly affect the shape of small voids. 

Figure [2] depicts the effect of cholesterol on (cos0}v 
vs. V. At all X, the voids become on the average slightly more 
perpendicularly oriented with an increasing V. An increas- 
ing cholesterol concentration makes voids of all sizes adopt 
average orientations closer to the direction of the bilayer nor- 
mal. We can therefore conclude that an increasing cholesterol 
concentration affects the orientation of at least the voids with 
4 X 10"^ nm^ <V < 0.13 nm^. 



H. Local Effects 

We have seen that a finite cholesterol concentration changes 
the properties of voids in a bilayer. It is natural to ask whether 
or not the changes are more pronounced in the close vicinity 
of cholesterol molecules. In principle, investigating the local- 
ity of the changes cholesterol brings about is simple: we just 
have to single out the voids whose CMs are located near the 
CM of a cholesterol molecule, and compare their properties to 
voids whose CMs are nowhere near cholesterols. There are, 
however, a few catches that should be avoided. 

First, studying arbitrarily large voids is meaningless. The 
CM of a very large void may or may not be close to the CM of 
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a cholesterol, even though some parts of the void might be in 
the vicinity of one or several cholesterols. In any event, such 
large clusters cannot very well be said to be localized. Hence, 
we have confined ourselves to voids that are smaller than a 
cholesterol molecule, i. e., voids with V < Vcho\- 

An additional problem is choosing a suitable criterion for 
the proximity of a void to a cholesterol molecule. The crite- 
rion should not be too strict, leading to automatic exclusion of 
all large and perpendicularly oriented voids from the class of 
voids close to cholesterol. This can be ascertained by moni- 
toring the principal radii of voids and choosing the criterion 
based on these. The choice can be facilitated by restricting 
the study to Region 3, where the voids will, by definition, be 
in the same region as the steroid ring structures. This allows 
us to consider distances that are projected to the plane of the 
bilayer. 

We have chosen to focus on concentration x = 20.3%, 
and the penetrant radius has been set to r = O.lOnm. After 
investigating the principal axes of voids in this case, we have 
decided to regard voids with CMs within 0.8 nm of a CM of 
a cholesterol molecule, making sure that both are in the same 
monolayer, as proximate to cholesterol. This decision is based 
on two observations: the radius corresponding to the thickest 
part of a cholesterol molecule in our study is approximately 
0.3 nm,^^ and the vast majority of voids in our cuiTent setup 
have principal radii smaller than 0.5 nm. 

The results are shown in Fig.^] Panel (a) features the size 
distributions computed for voids in the vicinity of cholesterols 
and for those far from cholesterols. To facilitate comparison, 
these distributions have been normalized by the average num- 
ber of voids close to cholesterols in Region 3 and by the aver- 
age number of voids far from cholesterol in Region 3. The re- 
sulting normalized distribution P{V) is such that it gives unity 
when integrated over Voids with radii smaller than 0.1 nm 
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FIG. 13: Effect of proximity of cholesterol molecules at x = 20.3 % 
with r = 0.10 nm. The data for voids close to cholesterols are shown 
in black and those for voids far from cholesterols in gray, (a) Nor- 
malized void size distributions, (c) Piaxja-i). (b) Orientational dis- 
tributions. 



are unaffected by the proximity of cholesterol. Larger voids, 
however, are affected: the vicinity of cholesterol reduces the 
probability of finding a void with V k, 14hoi/10 by a little 
less than a factor of ten. Further, near cholesterols there are 
no voids with 0.1 nm'' <V < V^hoi- 

Is the size distribution for voids far from cholesterol coin- 
cident with the size distribution for Region 3 of a pure DPPC 
bilayer? Considering voids with V < Vchoi and normalizing 
by the total number of voids in each case, we obtain distribu- 
tions that can be compared directly. It turns out that the two 
are nearly identical. Hence, further away from cholesterols 
the size distribution of voids is unaffected by the presence of 
cholesterol and practically identical to the coiTesponding dis- 
tribution in pure DPPC. 

The effect of the proximity of cholesterol molecules on void 
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shape is portrayed in panel (b) of Fig.^] It is evident that the 
shape of voids does not depend on the vicinity of cholesterols. 
Orientation of voids near cholesterols and far from them is 
portrayed in panel (c) of the same figure. The voids that are 
located close to cholesterols tend to be more oriented along 
the bilayer normal, while the ones far from cholesterols show 
a broader orientational distribution. 



V. CONCLUDING DISCUSSION 

We have studied the impact of cholesterol on the voids 
in simulated bilayers consisting of DPPC. Using data from 
atomistic molecular dynamics simulations, we have investi- 
gated the behavior of the total free volume with an increasing 
cholesterol concentration, extracted the void size distributions 
in different parts of the bilayer, and considered the effect of 
cholesterol on the shape and orientation of the voids. We have 
also probed whether or not the effects of cholesterol on the 
properties of voids are localized to the near vicinity of choles- 
terol molecules. In interpreting the data, we have, in the spirit 
of Ref. 9, made use of three-dimensional isotropic percolation 
theory. 

Our findings support the reduction of the total free volume 
with an increasing cholesterol content. The effects of choles- 
terol on individual voids appear to be the most prominent 
in those parts of the bilayer where cholesterol's steroid ring 
structures and the DPPC tails are located. Here cholesterol 
reduces the connectivity of free volume, thus altering the void 
size distributions. These changes are perceptible already at 
low cholesterol concentrations and become pronounced with 
approximately 20 % cholesterol in the bilayer. For instance 
a concentration of 29.7 % leads to complete removal of the 
voids that are approximately of the same size as a choles- 
terol molecule. An increasing cholesterol concentration also 
makes the voids in the region containing the steroid rings more 
elongated, and transforms the isotropic orientational distri- 
bution of voids present in pure DPPC to an anisotropic one 
favoring orientation of voids along the bilayer normal. All 
these effects appear to be localized to the proximity of choles- 
terol molecules, rather than being uniformly distributed to the 
whole bilayer 

Reasons for the changes in the region with the choles- 
terol ring structures and the hydrocarbon tails of DP- 
PCs can be found in previous computational studies on 
DPPC /cholesterol systems . 2 2.2^,49^50,5 1 Cholesterol increases 
the order parameters of the hydrocarbon tails by simultane- 
ously reducing their tilt and the number of gauche-defects 
present in them. The straightening of the tails should affect 
both the void size distributions and, in particular, the shape 
and orientation of the voids. This explanation is supported by 
the fact that the effects of cholesterol are confined to the prox- 
imity of cholesterol molecules: the effect of cholesterol on 
the ordering of the phospholipid tails is much stronger in the 
vicinity of cholesterol molecules. In addition to the order- 
ing, there is another, more obvious factor, which should cer- 
tainly influence the void size distributions. The bulky choles- 
terol molecule simply fills empty spaces of its own size in the 



bilayer, as suggested by Almeida et al. 

Our results for pure DPPC are in good qualitative agree- 
ment with the findings of Marrink et al.,^ although their study 
has been performed at T = 350 K and is based on 80 ps 
of data. The form of our void size distributions in different 
regions and for different permeants is in surprisingly good 
agreement with theirs. Their measure for the shape of the 
voids is different from ours, not enabling them to fully distin- 
guish between elongated and fractal voids. In addition, they 
study voids of all sizes, while we have restricted our studies on 
the shape and orientation of voids to the range 4 x 10~^ nm'^ 
<V < O.lSnm"^. Their qualitative conclusion is that larger 
voids may be more elongated and /or more fractal than small 
ones. We, on the other hand, find that the shape of the voids 
depends only weakly on the size, and that the voids become 
slightly less elongated with an increasing volume. As for ori- 
entation, despite the different measures, their conclusions are 
in accord with ours. They observe that the orientational dis- 
tribution of voids is most anisotropic is in Region 3, while 
Region 1 features an almost completely isotropic distribution. 

We also agree with Marrink et al.'' in that the complexity 
of the properties of the voids in a phospholipid bilayer im- 
ply that simple free area theories cannot be used for quanti- 
tative predictions. However, free volume arguments may be 
used for qualitatively explaining diffusion and permeation in 
bilayers. An increasing cholesterol concentration in DPPC bi- 
layers leads to monotonously decreasing lateral diffusion co- 
efficients for both DPPCs and cholesterolsi^^ Based on our 
current findings on the properties of voids, this is to be ex- 
pected. Intuitively, smaller and fewer voids, which are ori- 
ented along the bilayer normal, should slow down lateral dif- 
fusion of lipids and sterols. The details of what role the voids 
play for the various diffusion and permeation mechanisms in 
bilayers are challenging questions that large-scale simulation 
studies may elucidate in the near future. 

In addition to the static properties of voids considered in 
this study, also the dynamic properties should be important 
to dynamic processes in bilayers. It would be interesting to 
study the dynamics of void formation and the movement of 
individual voids. Such studies should give us further insight 
into mechanisms of lateral diffusion and solute permeation in 
membranes. 
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